% fig2c1_H.m ——— 绘制 Fig.2(c1)：<H> vs alpha
clear; clc; close all;

%% —— 模型及仿真参数 —— 
a     = 0.5;       % 论文中 a'
b     = 0.215;     % 论文中 b'
beta  = 0.2;       % 论文中 β
A     = 0.95;      % 强迫幅值
omega = 0.61;      % 固定 ω = 0.61

T     = 800;       % 总仿真时间
dt    = 0.01;      % 积分步长
T0    = 600;       % 丢弃前 transient 时间

%% —— 参数扫描 —— 
% 为了贴近论文截图，这里 α 从 1.6 到 3.2，共 200 点
alpha_list = linspace(1.6, 3.2, 200);
N = numel(alpha_list);
H_vals = zeros(1, N);

X0 = [0.01; 0.02; 0.02];           % 初始条件
tspan = 0:dt:T;

for k = 1:N
    alpha = alpha_list(k);
    % 1. RK4 积分
    X = rk4(@(t,X) rhs(t,X,a,b,beta,A,omega,alpha), tspan, X0);
    % 2. 丢弃 transient
    idx = tspan > T0;
    X_tr = X(:, idx);
    % 3. 计算能量时间序列
    H_series = compute_energy(X_tr, alpha, a, b);
    % 4. 取平均
    H_vals(k) = mean(H_series);
end

%% —— 绘图 —— 
figure;
plot(alpha_list, H_vals, '-','LineWidth',1);
xlabel('\alpha', 'FontSize',12);
ylabel('\langle H\rangle', 'FontSize',12);
title('平均能量 \langle H\rangle vs \alpha (\omega=0.61)', 'FontSize',14);
axis([min(alpha_list) max(alpha_list) min(H_vals)*0.9 max(H_vals)*1.1]);
grid on;

% compute_energy.m ——— 按论文式 (6) 计算能量 H
function H = compute_energy(X, alpha, a, b)
    % X: 3×M 矩阵，每列是一个时刻的 [x; y; z]
    x = X(1,:);
    y = X(2,:);
    z = X(3,:);
    H = 0.5*alpha.*x.^2 + 0.5*y.^2 + 0.5*a.*z.*x + 1.5*b.*z.^3.*x;
end
